log using "$Rep_smokelabor/2_analysis/output_logs/figure5.log", replace 

** panel a 
if 1 {
	
	use "$Rep_smokelabor/1_build/regdata/agegrp_county_quarter.dta", clear 
	
	foreach grp in A04 A05 A06 A07 A08 {
	preserve    
		
		keep if agegrp=="`grp'"
		tsset fe_countyqtroy rfrnc_yr
		
		* first diff: y(t) minus y(t-1)
		foreach v of varlist pc_qwi_payroll {
			gen d_`v'=`v'-L1.`v'
		}
		
		ivreghdfe d_pc_qwi_payroll (pm25=hms_deep) [aw=seer_pop] , a(fe_countyqtroy fe_styr) cluster(countyfip fe_stqtros)
	
	restore 	
	}
	
	** plot 
	import excel using "$Rep_smokelabor/2_analysis/output_figures/figure5.xlsx", sheet("age") clear first
	
	label define agegrp_lb 1 "age 25-34" 2 "age 35-44" 3 "age 45-54" 4 "age 55-64" 5 "age 65+", replace
	label values agegrp agegrp_lb
	
	gen pct_inc=b_inc*100/m_inc
	
	tw 	bar pct_inc agegrp, yaxis(1) col(teal%30) lcol(white) lw(0.001) barw(0.4) || ///
		rcap l_inc u_inc agegrp, col(gs9) lp(dash) yaxis(2) || ///
		scatter b_inc agegrp, col(blue) yaxis(2) ///
		ylab(-350(70)0, nogrid axis(2)) ylab(-3(0.6)0, axis(1) nogrid) ///
		xtitle("") ytitle("") ///
		yline(0, lcol(black) lp(dot) lw(0.6)) ///
		xlab(1 2 3 4 5,valuelabel labsize(small)) ///
		legend(col(2) size(small) order(1 "% change (left axis)" 3 "\$ change (right axis)")) ///
		graphregion(col(white))
	gr export "$Rep_smokelabor/2_analysis/output_figures/figure5_a.pdf", replace
								
}

** panel b
if 1 {
	
	use "$Rep_smokelabor/1_build/regdata/naics2_county_quarter.dta", clear 
	
	forv grp =1/20 {
	preserve    
		
		keep if g_industry==`grp'
		tsset fe_countyqtroy rfrnc_yr
		
		* 1-year change: y(t) minus y(t-1)
		foreach v of varlist pc_qwi_payroll {
			gen d_`v'=`v'-L1.`v'
		}
		
		ivreghdfe d_pc_qwi_payroll (pm25=hms_deep) [aw=seer_pop] , a(fe_countyqtroy fe_styr) cluster(countyfip fe_stqtros)
	
	restore 	
	}
	
	egen fe_indcountyqtroy=group(industry countyfip rfrnc_qtroy)
	tsset fe_indcountyqtroy rfrnc_yr
	
	* first diff: y(t) minus y(t-1)
	foreach v of varlist pc_qwi_payroll {
		gen d_`v'=`v'-L1.`v'
	}
	
	compress g_industry
	local cmd ivreghdfe OUTCOMEVAR  (pm25=hms_deep) [aw=seer_pop] , a(fe_countyqtroy fe_styr) cluster(countyfip)
	wyoung d_pc_qwi_payroll, cmd(`cmd') familyp(pm25) bootstraps(100) subgroup(g_industry) seed(20) cluster(countyfip) replace
	saveold "$Rep_smokelabor/2_analysis/output_figures/figure5_b_westfall_young.dta", replace 
	
	** plot 
	import excel using "$Rep_smokelabor/2_analysis/output_figures/figure5.xlsx", sheet("ind") clear first

	label define industrygrp_lb 1 "Agriculture" 2 "Mining" 3 "Utility" 4 "Construction" 5 "Manufacturing" 6 "Wholesale" 7 "Retail" 8 "Transportation" 9 "Information" 10 "Finance" 11 "Real estate" 12 "Professional" 13 "Management" 14 "Admin" 15 "Education" 16 "Healthcare" 17 "Entertainment" 18 "Accommodation" 19 "Other services" 20 "Public admin" , replace
	label values industrygrp industrygrp_lb

	gen pct_inc=b_inc*100/m_inc
	gsort pct_inc
	gen s_industrygrp=_n
	
	tw 	bar pct_inc s_industrygrp,  horizontal yaxis(1) col(teal%30) lcol(white) lw(0.001) barw(0.4) || ///
		rcap l_inc u_inc s_industrygrp, horizontal col(gs9) lp(dash) xaxis(2)   || ///
		dot b_inc s_industrygrp, dotex(no) ndots(0) horizontal col(blue) xaxis(2) || ///
		dot b_inc s_industrygrp if inlist(industrygrp,1,2,9,10,12,13,15,17), dotex(no) ndots(0) horizontal col(white) msize(0.8) xaxis(2)  ///
		xline(0, lcol(teal*0.5)) ///
		ytitle("") ///
		xtitle("Coefficient") ///
		legend(col(2) size(small) keygap(0.5) region(lcol(white)) order(1 "% change (bottom axis)" 3 "# change (top axis)")) ///
		ylab(1(1)20, nogrid valuelabel angle(0) labs(vsmall)) ///
		xlab(-4(1)1, axis(1)) xlab(-24(6)6, axis(2)) ///
		xtitle("", axis(1)) xtitle("", axis(2)) ///
		graphregion(color(white)) ///
		ysize(5.5) yscale(reverse)
	gr export "$Rep_smokelabor/2_analysis/output_figures/figure5_b.pdf", replace
}

log close
